基于核平滑构建物种分布范围：

代码来源

## 基于核平滑计算物种分布范围：
https://jamesepaterson.github.io/jamespatersonblog/04_trackingworkshop_kernels

案例：

##基于喜旱莲子草的物种分布数据构建内核分布：
as <- read.csv('D:/xh2/f1_points/xhas.csv')[,2:3] %>% na.omit() %>% round(4)
tail(as)
library(sp)
plot(as2)
map(add =T)
as2 = SpatialPoints(cbind(as$longitude, as$latitude), proj4string = CRS("+init=epsg:4326" ))
as.sp <- spTransform(as2, CRS("+init=epsg:3857"))

as.sp@coords
## 基于登高线构建核密度分布：
library(adehabitatHR)

kernel.ref <- kernelUD(as.sp, h = "href", grid = 1000,hlim = c(0.5, 1.5))  # href = the reference bandwidth
image(kernel.ref) # plot
## 将等高线栅格化：
kde <- raster(kernel.ref)
plot(kde)
# 可视化：
library(tmap)
tm_shape(kde) + tm_raster("ud")

## 构建等高线：
as.kernel.poly <- getverticeshr(kernel.ref, percent = 90) 
plot(as.kernel.poly)
points(as.sp)


##  将构建的范围转为wgs84格式的shp并导出：
zk_area<- spTransform(as.kernel.poly, CRS("+init=epsg:4326"))
library(rgdal)
writeOGR(zk_area, "C:/Users/admin/Desktop/aa.shp",layer="ll" , driver="ESRI Shapefile")

## 读取shp文件：
library(rgdal)
rgdaltest<-readOGR('C:/Users/admin/Desktop/l2.shp')
rgdaltest@proj4string
plot(rgdaltest)
coordinates(as) <- c("longitude","latitude" )
proj4string(as) <- CRS("+init=epsg:4326")
points(as)


## 基于范围边界的核密度图估计(非截取)####
library(sp) # Spatial objects
library(rgdal) # If you want to read boundary as a shapefile
library(adehabitatHR) # For kernels
library(dplyr) # For filtering data
library(magrittr) # For piping
library(xlsx)
st1 <- read.xlsx("mtcars.xlsx", 1)


turtles <- read.xlsx(file = "C:/Users/admin/Desktop/tracking_sample.xlsx",
                    1) %>%
  filter(id %in% c("T003", "T004"), # we will only use 2 individuals in example (T003 and T004)
         !is.na(x), # Remove NA coordinates
         !is.na(y)) 

## 基于阻碍构建数据核密度分布：
# 需要注意这和截断连续密度分布不同；
# 其中boundary必须使用线文件；

bound.kernel.t003 <- turtles.sp[turtles.sp$id == "T003","id"] %>%
  kernelUD(.,
           # h must be numeric 
           # here I give the numeric argument for "href" from our estimate above
           h = trad.kernel$T003@h$h, 
           # If you don't set grid, the default often won't follow the boundary
           # I set grid to have 1000 x 1000 pixels
           # You may have to adjust (will depend on your scale)
           grid = 1000, 
           # The boundary argument takes SpatialLines objects
           boundary = bound) 

# Get the vertices for T003 and extract as a polygon
bound.kernel.t003.poly <- bound.kernel.t003 %>%
  getverticeshr(., percent = 95)

# Make a kernel with boundary for T004
bound.kernel.t004 <- turtles.sp[turtles.sp$id == "T004","id"] %>%
  kernelUD(.,
           h = trad.kernel$T004@h$h,
           grid = 1000,
           boundary = bound) 

# Get vertices for T004
bound.kernel.t004.poly <- bound.kernel.t004 %>%
  getverticeshr(., percent = 95)